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ABSTRACT 

We use numerical simulations to investigate the cusp at the centre of elliptical 
galaxies, due to the slow growth of a super-massive black hole. We study this problem 
for axisymmetric models of galaxies, with or without rotation. The numerical simula- 
tions are based on the 'Perturbation Particles' method, and use GRAPEs to compute 
the force due to the cusp. We study how the density cusp is affected by the initial 
flattening of the model, as well as the role played by initial rotation. The logarithmic 
slope of the density cusp is found to be very much insensitive to flattening; as a conse- 
quence, we deduce that tangential velocity anisotropy- which supports the flattening- 
is also of little influence on the final cusp. We investigate via two different kinds of 
rotating models the efficiency with which a rotation velocity component builds within 
the cusp. A cusp in rotation develops only for models where a net rotation compo- 
nent is initially present at high energy levels. The eventual observation of a central 
rotational velocity peak in E galaxies has therefore some implications for the galaxy 
dynamical history. 
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1 INTRODUCTION 

HST observations have shown that a central density 
cusp is present in most, if not all, elliptical galaxies (Lauer 
et al. 1995, and companion papers). Furthermore the slope 
of the density cusp has shown a dichotomy between high 
mass (luminosity) systems and low mass (luminosity) ones, 
analogous to what has been found for other observational 
properties (Bender et al. 1989). Low mass ellipticals tend to 
host steep cusps, the radial profile of their luminosity p(r) 
having a logarithmic slope 7 = —dlog(p) /dlog(r) ~ 1.9±0.5. 
On the other hand, luminous ellipticals have shallower cusps, 
with 7 ~ 0.8 ± 0.5 (Gebhard et al. 1996). 

The most popular explanation for the formation of a 
density cusp is a central super-massive black hole (BH)- see 
Richstone et al. (1998), Kormendy & Richstone (1995) for 
reviews. Such BHs are believed to be present in a large frac- 
tion (possibly close to 1, Haenelt & Rees 1993, Tremaine 
1997) of present day galaxies. Their influence on a col- 
lisionless galactic nucleus was first addressed by Peebles 
(1972), Young (1980), Goodman & Binney (1984) using 
semi- analytic models. 

All these models suppose an isolated galaxy, at the cen- 
ter of which a BH grows by gas accretion. The slope of the 



cusp produced in such models is within the range observed 
for the less massive E galaxies. In this paper we shall also 
consider the formation of a cusp in an isolated galaxy, as we 
are interested in the origin for the steep cusp (7 ~ 1.9) of 
low luminosity ellipticals. On the other hand, numerical sim- 
ulations (Makino & Ebisuzaki, 1996, Quinlan & Hernquist 
1997, Nakano & Makino 1999) have shown that a BH sink- 
ing within the core of an elliptical galaxy, or a binary of BHs 
produce a shallow cusp, similar to that observed for massive 
E galaxies. Such a BH binary could result from a merger. 
This picture is therefore compatible with the general belief 
(e.g. Nieto & Bender 1989) that the apparent existence of 
two classes of ellipticals is related to a more important role 
played by merging and interactions in the history of massive 
E galaxies. 

All the semi-analytic models of cusp formation due to a 
central BH use the adiabatic invariance of actions to derive 
the distribution function after a BH has slowly grown (see 
Young 1980). Since actions are explicitly known for spheri- 
cal systems, analytical models have been derived in spher- 
ical symmetry. Based on such quasi-analytic computations 
for spherical models, Quinlan et al. (1995) investigated the 
consequences of the initial density profile on the slope of the 
cusp, and addressed the influence of a velocity anisotropy. 
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For non-spherical systems - except for Stackel poten- 
tials - an expression for the actions is unknown, and there- 
fore models based on the conservation of actions have not 
been computed hitherto. In this paper, we propose to use 
numerical simulations with high central resolution, based 
on the PP method (Leeuwin et al. 1993), to investigate the 
cusp formed in axisymmetric systems by the growth of a 
BH. 

Many low luminosity galaxies display little evidence for 
triaxiality, but are close to axisymmetry; they can be mod- 
eled as two-integral models, with f(E,L z ). Such a distri- 
bution function (d.f.) implies velocity dispersions obeying 
or = °"z, meaning that flattening would be due to an ex- 
cess of azimuthal motion (which increases the net rotation, 
or the tangential anisotropy). This does not preclude that 
three-integral models would not do as well, or better (e.g. 
review by Merritt, 1998). Many of those galaxies however 
exhibit kinematic features consistent with isotropic rotators 
models: at least for them, an excess of tangential motion is 
plausibly the main support for their flattening. 

The BH itself may have perturbed orbits sufficiently to 
erase the part of memory for initial conditions that corre- 
sponds to conservation of a third integral (Norman et al. 
1985, Gerhard & Binney 1985, Merritt & Quinlan 1998). In 
this case axisymmetry is a consequence of the BH growth. 
Nevertheless, one still ought to investigate the scenario in 
which, for some of these galaxies, the presently observed ax- 
isymmetry already existed before the BH growth. 

This is the goal of this paper, where we will study the 
cusp generated by the growth of a central BH within ax- 
isymmetric systems having various degrees of rotation. If a 
BH grows in a roughly axisymmetric galaxy, do the initial 
flattening and rotation have a sizeable effect on the evolu- 
tion, so that it deviates from that of the well known spheri- 
cal case? Could we in certain cases be able, eventually with 
more detailed observational data, to infer from the proper- 
ties of presently observed cusps any information about the 
pre-black hole galaxy? Should one be cautious in certain 
cases about deriving the BH mass using the spherical adia- 
batic model? 

This paper is organized as follows: first, the existing 
models are briefly recalled in section 2. We then specify the 
initial conditions (section 3) and the numerical techniques 
used for this work (section 4). Results for non-rotating and 
rotating models are given in sections 5 and 6, respectively. 
We conclude in section 7. 



2 PREVIOUS MODELS FOR GALACTIC 
CUSPS 

Previous models are based on a few basic assumptions, and 
ours will follow suit. The growth of the BH is supposed to be 
due to gas accretion, and thus does not deplete the original 
stellar component. The existence of the gaseous component 
is otherwise neglected, so that the models follow the evolu- 
tion of a purely stellar component. The growth time tsu of 
the black hole is taken to be at maximum ~ 10 7 Yrs. This 
is long enough with respect to the stellar dynamical times, 
so that we can compare our models with the adiabatic com- 
putations, but not too long for obvious computation time 
reasons. It is furthermore reasonnable if compared to the 



'Salpeter time', or timescale for the mass accretion onto a 
black hole: 

t s = Mbh/Mbh ^ e x 4 x 10 8 Yrs. (1) 

In the above, e is the radiative efficiency, which is most often 
e < 0.1 in current models of accretion disks (cf. Richstone 
et al. 1998), or in estimates derived from quasars number 
counts (van der Marel 1997). 

Furthermore the models neglect collisions. It has been 
estimated in the past that only some of the spiral galaxies 
bulges are dense enough to experience significant two-body 
interactions (see e.g. Kormendy 1988), although these esti- 
mates were based on the belief that elliptical galaxies had 
cores. Let us suppose for instance that up to ~ 0.1 — 1% of 
an E galaxy mass, corresponding to say TV ~ 10 6 — 10 9 stars, 
is enclosed in the central ~ 10 pc, where velocity is of order 
a few times ~ 100 km/s. The crossing time in these regions 

should be T cr ~ 10 5 (^j^) ( ^j^/s) ) Yrs - Therefore the 

relaxation time would be (e.g. Binney & Tremaine 1987): 

t rel * N/(8lnN)t cr * 10 9 (^) (^) (&) Yrs. 

Stellar collisions are liable to play some role over the age of 
the universe in such a system with an already highly concen- 
trated density. For the models we consider here, however, a 
collisionless model is justified, given that the initial galaxy 
has a large core, and that the BH grows in ~ 10 7 Yrs. Nev- 
ertheless it should be kept in mind that collisional processes 
could play a role in the evolution of the center-most part 
after the cusp has developed. 

The spatial extent of the cusp is given approximately 
by the influence radius ri of the BH, which is the distance 
up to which the BH potential dominates. This is estimated 
as: 



where ao is the central velocity dispersion, and Mbh is the 
BH mass. 

2.1 Analytical models: spherical adiabatic models 

Analytical models assume that the central BH grows over 
time scales that are long compared to dynamical times in the 
central region. The final d.f. can then be deduced from adi- 
abatic conservation of actions (Peebles 1972, Young 1980). 

In an isotropic system having initially a core, the cusp 
is predicted to have a slope 7 = 3/2 (Young 1980; see also 
Tremaine 1997). This value, however, holds if the BH mass is 
less than roughly the core mass; when larger, the self-gravity 
of the cusp becomes important and the slope increases up to 
values ~ 2 (e.g. Young 1980). For centres that are initially 
'non-analytical' - i.e. whose density profile has a non-zero 
derivative at the origin- Quinlan et al. (1995) derive larger 
slopes 7, with higher values for stronger initial cusps. 

Goodman & Binney (1984) provide the following simple 
picture for the way orbits are transformed as the cusp grows. 
The fact that radial action is conserved plays an especially 
important role in the shape evolution of elongated orbits. 
It means that the radial excursion r max — r m i n does not 
change much, but the potential center, which is originally 
at the orbit center (roughly harmonic potential), is finally 
at one of its foci (keplerian potential). Thus for elongated 
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orbits eccentricity has decreased, which means also that the 
motion is more circularly biased. The adiabatic analytical 
models already cited predict indeed a moderate tangential 
anisotropy. This trend for rounder orbits is also present in 
numerical results by Norman et al. (1985), who studied how 
the shape of orbits close to a growing point mass is modified. 

Initial anisotropies in the stellar velocities appear to 
have a weaker influence than an initial cusp (Quinlan et al. 
1995). Due, however, to the small number of available ana- 
lytical anisotropic models, this problem could not be thor- 
oughly examined. One of the few anisotropic cases for which 
the cusp slope has been semi-analytically evaluated is a 
spherical model fully made of circular orbits. From conserva- 
tion of energy and angular momentum, the final logarithmic 
slope is found to be 7 = 9/4 (Quinlan et al. 1995 |j), signif- 
icantly higher than for the isotropic case (7 = 3/2). There- 
fore we could expect that in a model with some intermediate 
degree of tangentially biased motion, 3/2 < 7 < 9/4. Since 
our axisymmetric models are flattened in their centre by an 
excess of tangential motion, it is interesting to ask whether 
such an anisotropy will have a sizeable influence on the cusp 
slope. We may also ask whether a systematic rotation will 
bring any changes. These are some of the questions we will 
try to answer in this paper. 

For the very few galaxies which had been observed 
ten years ago with sufficient central resolution, observations 
showed rotation velocities rising towards the center. This 
prompted attempts to compute the amount of angular mo- 
mentum brought to the center by the forming cusp (Lee 
& Goodman 1989, hereafter LG89 ; see also Cipollina & 
Bertin 1994, hereafter CB94). These works, however, had 
to assume a spherical potential in order to derive analytical 
solutions. One aim of this paper is to address this problem 
anew without any such approximation. In sections 5 & 6, we 
investigate the influence of initial rotation for self-consistent 
axisymmetric models. Our motivations do not follow those 
of LG89, because higher quality, more recent observations 
no longer show evidence for more rotation than dispersion 
in the centre (see e.g. van der Marel et al. 1997, Joseph et 
al. 2000 for HST observations of M32, and Kormendy et al. 
1998, regarding NGC 3377). 



2.2 Numerical methods 

The pioneering simulation by Norman et al. (1985) did not 
search for any detailed feature in the center, but was meant 
instead to measure the shape changes caused at moder- 
ate and large radii by a central growing BH. Due to the 
limited numerical capacity available at the time, their re- 
sults displayed substantial particle noise, but showed al- 
ready the trend towards axisymmetrization. Merritt & Quin- 
lan (1998) recompute this problem with a considerably im- 
proved accuracy, obtained by using a larger number of par- 
ticles (N ~ 1.1 x 10 5 instead of N = 2 x 10 4 for Norman 
et al. 1985), and by enhancing the central particle density. 
To this end, they replace each strongly bound particle by 



* using their notations in which initially p(r) ~ r c , the value 

for the slope given in their Eq. (17) should actually be C = 3 — 

3 
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a number of less massive particles distributed over its orbit 
(see §1.4). 



Sigurdsson et al. (1995) built a numerical code to follow 
the growth of the cusp itself, and tested it for the spheri- 
cal case. Their code makes use of the self-consistent field 
method (Hernquist & Ostriker 1992). This is formulated as 
a perturbation method, in the sense that it uses a family 
of density-potential functions, of which the zero-order term 
corresponds, or is close, to the initial model. At each time 
step during the simulation, the coefficients of the expansion 
for the density yield the new potential. The less the system 
evolves, the fewer are the terms required in the expansion 
to provide a given level of accuracy. In the implementation 
by Sigurdsson et al. (1995), the number of particles in the 
central parts is enhanced by using a mass spectrum for the 
particles. This is achieved by considering for the particles 
distribution, instead of the initial d.f. fo(E), a distribution 
fo(E) /m(E, L), such that the mass m of a particle is less 
than unity if the pericenter of its orbit is smaller than some 
radius, and is unity if it is larger. This way, orbits reach- 
ing close to the centre are represented by more numerous, 
lighter particles. A numerical simulation based on this code, 
and using 512, 000 particles, has been performed by van der 
Marel et al. (1997), in order to check the stability of their 
dynamical model for M32, i.e. an axisymmetric model in- 
cluding a central supermassive BH. 

In the present paper, we use a numerical perturbation 
technique which allows to concentrate particles at the cen- 
ter of the model, without entailing at the same time high 
particle-noise for the potential at larger radii (§4). In addi- 
tion, we perform a subdivision of the central particles when 
necessary, to further increase the resolution at the center. 
These numerical simulations are used here in order to extend 
the models of BH induced cusps to non-spherical systems. 



3 THE INITIAL AXISYMMETRIC MODEL 

3.1 Axisymmetric, non-rotating model 

Fully analytical axisymmetric models for elliptical galaxies 
are scarce (Lynden-Bell 1962, Dejonghe 1986, Evans 1994, 
cf. discussion in Hunter & Qian 1993). Most of them are 
sums -often infinite sums - of Fricke terms. We take here 
the simple model by Lynden-Bell (1962), which uses only 
two Fricke terms: 



f (E,L z ) = FiE 7/2 + F 2 L 2 z E l:i/2 if£ > 



(3) 







otherwise 



where E is the binding energy E = <f>o — v 2 /2 and L z is 
the component of angular momentum about the symmetry 
axis. This distribution corresponds to a flattened Plummer 
potential: 



MR,*) = 



M 



{{R 2 + z 2 +a 2 p ) 2 -2b 2 p R 2 y/ A '' 



(4) 



where we have used cylindrical coordinates. In the above, Mq 
is the total mass of the model, ap and bp are two parameters, 
and Fi and F2 are normalizing numerical constants, which 
expression is given by Lynden-Bell (1962): 
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Figure 1. Left panel: isodensity contours of the Lynden-Bell 
model with ap = 4 and ftp = 2.2; the axis ratio for the density is 
shown as a function of the cylindrical radius in the insert. Right 
panel: velocity dispersion profiles cr|j = a\ and it? along the R 
axis for the same model 
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We take in the following ap — 4 and Mo = 1 in the units 
of the simulation. These are such that the mass unit is 
10 11 Mq, the length unit is 1 kpc, and G = 1. 

The core radius is usually defined as the radius at which 
the surface brightness £(i?, z) has decreased to half its cen- 
tral value So. This corresponds for the present model to 
R1/2 — 0.76 ap. The mass enclosed within the E = est 
isodensity level passing at (R1/2, z = 0) is Moore — 0.2 Mo. 

The flattening decreases with radius, so that the mod- 
els tend towards spherical symmetry at large radii. The 
larger bp, the larger the central flattening. b p however must 
be smaller than bp max — 0.64 a p otherwise the isodensity 
contours have a dimple on the z-axis. We shall consider a 
model close to the model with maximum flattening, by tak- 
ing bp — 0.55 ap. For this model the axis ratio is about 
constant, with value c/a ~ 0.55 (Fig. [l]) within the region 
where the cusp will develop (see above the definition of the 
influence radius), 

Since the d.f. depends only on E and L z , the velocity 
dispersions obey ap — a z . The flattening is sustained by an 
excess of motion in the <f> direction, with respect to motion 
in the Rot z directions. If velocity fields that are moderately 
biased towards tangential motion can affect the final density 
cusp, we should witness it in our computations. 

The distribution function in Eq. H is even in L z and cor- 
responds to a non-rotating system. In this case (c/. Lynden- 
Bell 1962): 



2 
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(6) 



(7) 



with C = 5bp(2ap - bp)/(47rM s ). 

Fig. [l] displays the isocontours for the density, as well 
as the profile along R for the velocity dispersions along each 
direction. 



3.2 Rotating models 

Models with positive net rotation are obtained by changing 
the sign of L z for a fraction of the orbits having negative 
L z . If 77 is the fraction of orbits being reversed, we obtain: 



f(E,L z ) = 



f (E,L z )(l- V ), 
fo(E,L z )(l+r,), 



L z 



< 
> 



In order however to avoid a discontinuity at L z 
distribution function must have f(E,L z ) = f{E, 
L z — > 0. We therefore take for initial model: 

f5(E,L,) = f (E,L,) x {1 + tw(L,)), 



(8) 

the 
! for 

(9) 



where p{L z ) is some function having value ±1, except for 
a small range around the origin, where it goes to zero. A 
simple choice is: 



p{L z ) = 



sin te) 

sign(L z ), 



\L Z I < Lrn 

otherwise 



(10) 



Here L m is a parameter which determines the extent of the 
L z region where rotation is smaller than r\. Its choice is of 
special interest for what follows. 

The simplest choice for L m is a constant, independent 
of energy (c/. CB94). However, by taking some function 
Lm(E), we can build models where, for any given energy, 
a constant fraction of orbits contribute to rotation. This is 
similar to the case studied by LG89. We investigate both 
types of models, as we know from comparing these two pa- 
pers that they yield different results. 

We thus consider two kinds of initial distributions, that 
use respectively: 



Pi(L z ) 
P2(L Z ,E) 



p(L z ) withL„ 
p(L z ;E) withL„ 



xL c (E), < x < 1 (11) 



In the above definition, L C (E) is the angular momentum 
for the circular orbit with energy E, whereas Lq is some 
cons tant, w hich is conveniently expressed in units of lo = 
ap\/*o(0). 

The mean tangential velocity is given by: 

/ d 3 vfS(E,L z )(L z /R) 



v (f> (R,z) = 



Ve = 



J^vfS(E,L z ) 
The effective rotation is measured by: 
/ d 3 rd 3 vf(E,L z )L z 
J d 3 rd 3 vf(E,L z )\Lz\ ' 



(12) 



(13) 



which is in practice slightly smaller than 77 in the definition 
of /o given by Eq. M (and exactly 77 in the limit L m —> 0). 



4 NUMERICAL METHOD 

A crucial feature for any A^-body code aimed at studying the 
density cusp is of course the central resolution. Most parti- 
cles should ideally orbit within a few times the 'influence 
radius' rj of the BH (Eq. ^). In this paper, we use a scheme 
which is formulated as a perturbation method, and allows 
to choose freely the sampling in phase-space. The main ad- 
vantage of our scheme is that even if we have few particles 
outside n, where evolution is weak, the potential is fairly 
well represented (since it is dominated by the analytical $0 
term, which has zero particle-noise). 
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4.1 The 'perturbation particles' scheme 

The PP scheme has been described in detail elsewhere 
(Leeuwin et al. 1993). In this scheme, one writes the colli- 
sionless Boltzmann equation (CBE) as a perturbation prob- 
lem. The unperturbed state corresponds to the initial an- 
alytical model, while the perturbation corresponds to the 
evolution of the system. Thus, the system d.f. is formally 
written as / = fo + fi{t), where fi(t) is found by integrat- 
ing explicitly the CBE along orbits in an N-body realization. 
The choice of the orbits to be integrated during the run - 
that is to say the initial distribution of particles /s(r, v)- 
is in principle arbitrary, but one should aim at a good sam- 
pling of phase-space regions where large perturbations are 
expected. The function fs is chosen to be a function of the 
system classical integrals (here E,L Z ), so that in the absence 
of time evolution the sampling of phase-space is stationary. 

If the system potential varies with time, fs can not 
be made a stationary function. However, both fs and / 
obey the CBE equation. Therefore, noting the phase-space 
coordinate w = (r,v), we may simply record the initial 
value fs{w(to)) for each particle, since along each orbit 
fs{t) = fs(to)- The mass attributed to the particle running 
on that orbit is weighted according to the local phase-space 
density: 



mi(t) 



/sMto))' 



(14) 



as explained in Leeuwin et al. (1993) and Wachlin et al. 
(1993). 

On the other hand, the CBE is explicitly solved for f\: 
the accuracy of orbits integration can then be monitored 
through the relative variations of /. The perturbation along 
an orbit is given by: 



d/i(t) 
At 



[fo, <&BH + $s 



(15) 



where &bh = Mbh/t, and $ s . 9 . is the self- gravitating po- 
tential of the cusp. 

The computation of moments for /i amounts in fact to 
a Monte-Carlo integral with sampling fs- For instance, the 
perturbed total mass is evaluated as: 



Mi = / fi{w) a w ~ S i=1 - 



fs(m) 



L i=1 mii 



(16) 



To compute local values of the density field, or of any other 
field, the sum over the particles is of course limited to a small 
region of space where the quantity of interest is roughly 
constant. Moments of the total distribution / = fo + fi are 
given by: 



< V" >= 



fifo + h) V" d 6 w 
/(/o + /i)d 6 ™ ' 



(17) 



where V represents any component of the velocity vector. 
This is estimated as: 



<V n > (r) 



/ fo V" cfw + //if d 6 



po(r) +pi(r) 

,h 



po <v n > + x t 4±Mv? 



po < V n >o + E»mnV7 
po + pi 



(18) 



The sampling distribution fs is most efficient for the 
evaluation of the above Monte-Carlo integrals, when it is 
most similar to the integrants. The integration accuracy can 
be checked by monitoring the total mass in the perturbation, 
which should remain equal to zero. A convenient measure for 
the associated error is Mi/||Mi||, where we note ||Mi|| = 
Ei|mii|. 

The self-gravity of the cusp corresponds to the force due 
to the perturbation mass. Thus, a particle j will be subject 
to a self-consistent force term: 



r; - r,- 



(19) 



po + pi 



4.2 Orbit integration 

Close to the BH, the force gradients are of course very large. 
Some orbits therefore demand extremely short time-steps 
during small time intervals. For such problems, adaptative 
time-steps is the most convenient scheme. In a ID simulation 
(Leeuwin & Dejonghe 1998) that we ran to check whether 
the PP method could be advantageous for the present prob- 
lem, we have experimented block time-steps integration (us- 
ing either a leap-frog, or a RK4 scheme). In a block time-step 
scheme, the particles are sorted into groups at regular time 
intervals, and each group advanced with its smallest time- 
step until the next sorting out. This makes the scheme a 
costly one, since a set of particles may be advanced with a 
needlessly tiny timestep. 

Symplectic methods (Wisdom & Holman 1991, Saha & 
Tremaine 1992) , because they conserve phase-space volume, 
would in principle be useful to guarantee conservation of our 
sampling phase-space density (Binney, private communica- 
tion). Unfortunately, symplectic integrators with adaptive 
time-steps do not yet exist (see Duncan, Levison & Lee 1998 
for a version with block time-steps). 

We use in this work the ODEINT routine due to Press 
et al. (1986), which is an adaptative individual time-steps 
scheme, based on a 4th-order Runge-Kutta time interpola- 
tion. The conservation of / along orbits is better than 10 -4 
over each entire simulation. 



4.3 Force computation with GRAPE machines 

The BH mass was updated at each time-step. This is 
straightforward since the time dependence is explicit. More- 
over, within the BH influence radius where the BH force 
term dominates, a precise evaluation is required. 

The self-gravitating force (Eq. [llj) is computed by direct 
summation of the individual particles contributions, using 
a GRAPE machine. A small modification of the standard 
software was necessary since the PP masses can be either 
positive or negative, a possibility not foreseen for GRAPE . 
Furthermore, timesteps within the cusp are very small, so 
that the cost of re-evaluating by direct summation the inter- 
particles force at each time-step for the most bound orbits 
would be prohibitive. Since the BH grows very slowly (~ 10 7 
Yrs) com pared to the central dynamical times (~ 10 5 Yrs, 
see §2.1), the force field due to the cusp also changes slowly 
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with respect to individual crossing times. Moreover, it is neg- 
ligible at large radii where dynamical times are larger (see 
Fig. |^). Therefore the self-consistent force field was eval- 
uated using GRAPE only at regular time intervals At. To 
derive the self-gravity force experienced by each particle at 
intermediate times, we proceed as follows: 

• The self-consistent forces are evaluated using GRAPE 
at regular time intervals At. This yields the three compo- 
nents (F X i, F y i, F z i) of the force for each particle i. We infer 
for each particle the two components (Fm,F z i) along R and 
z in cylindrical coordinates. 

• Using a cell-in-cloud (CIC, cf. Hockney & Eastwood 
1981) scheme in polar coordinates, we compute the force 
field (Fr,F z ) on a grid in R, z (we average over the az- 
imuthal angle (f>). 

• At each time-step within At, the force (Fr, F z )(ri(t)) 
at the current position of each particle is interpolated from 
the grid values (with again a CIC). 

The grid has 20 x 20 cells, and extends in both directions 
from x\ — rj/200 to X20 = 4a p . Points are spaced with a 
logarithmic increment: 



Xi+l 



%i -I- dXi k , 



(20) 



where a and k are two real parameters numerically deter- 
mined after the choice of x\,dx\ and £20- We have taken 
dxi — x\. 

For the force evaluation with GRAPE we use a softening 
length e = 10~ 4 ap. The force is computed 1000 times during 
the BH growth. At each time-step, the value for the self- 
gravity force is linearly interpolated from the grid values, 
using again a CIC scheme. 

A run using ~ 3 x 10 J particles, and the parameters 
for the standard case defined hereafter (§5.1), takes approx- 
imately 2 days on the Marseille 5 board GRAPE-3AF sys- 
tem, and somewhat longer on the GRAPE-4 system with 41 
chips (see Athanassoula et al. 1998 for the characteristics of 
the Marseille GRAPE -3AF system, and Makino et al. 1997, 
Kawai et al. 1997 for the characteristics of the GRAPE -4 
system). For = 300, 000 particles, one force evaluation by 
GRAPE -3AF takes approximately 140 s, therefore ~ 38h in 
the simulation are devoted to the computation of the self- 
gravity force field. 

4.4 Choice of the sampling distribution 

In order to ensure a stationary distribution (in the absence 
of any perturbation), f s must be a function of the isolating 
integrals of motion for the unperturbed potential $o- We 
are limited in the choice of analytical distribution functions 
f s {E, L z ), as we were for fo(E, L z ), because such models are 
scarce in the iterature. The known d.f.'s for axisymmetric 
systems are in general expressible as Fricke series, for which 
the self-consistent density can be calculated. The simplest 
choice is thus to consider a term from a Fricke series 



f s (E, L z ) = FE a L 



a T 



(21) 



where E = $0 — \ and F is a normalizing constant. The 
corresponding density (see e.g. Dejonghe & de Zeeuw 1988) 
is given by: 

J/3 ^4-/3/2+3/2 



p s (R,z) = AF R $° 



(22) 
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Figure 2. Histogram for the normalized distribution of particles 
for each of the two sampling distributions fgi and f^2 discussed 
in the text, as a function of the cylindrical radius R (left) and of 
the spherical radius r (right) 



with A = r ^p ( ^ 1 _^y^ ) 2 1 ^2) 2 ' > • We can not generate a density of 
particles with a strong cusp using such a term. Models with a 
cusp in R can in principle be built by choosing — 1 < f3 < 0, 
however (i) their cusp in 71 has a rather faint logarithmic 
slope —1 < f3 < 0; (ii) these are d.f.'s with a cusp in L z , so 
their sampling of velocity space is very inhomogeneous. This 
is actually less unfavourable than it seems, since it favours 
nearly-radial orbits, which are actually th ose which expe- 
rience the strongest perturbation (see §5.5). More freedom 
is available regarding the density decay at large radii. A 
steeper decline of density, with respect to the unperturbed 
model, follows from taking a > 7/2 + 13/2, as can be easily 
verified. We are not worried by the low corresponding par- 
ticle density at large radii, since the potential there will be 
dominated by the analytical term. 

Here we experiment with two very different Fricke terms 
for f s , in order to make sure that the choice of the sampling 
does not affect the measure of the cusp slope, which is our 
main objective. We use: 



(i) /si: p = 0, a = 8. 

This is a model with a core, but its density falls much faster 
than the mass density of the unperturbed model. The central 
particle density is enhanced by splitting particles according 
to their energy and angular momentum, in a way similar to 
Merritt & Quinlan (1998). First particles are binned accord- 
ing to their binding energy Eq. For the most bound particles, 
we evaluate their pericenter Rmin{E, L z ). If R m in is smaller 
than the final BH influence radius (see §4.3), the orbit is se- 
lected for splitting. It is then integrated, within the unper- 
turbed potential, for sufficient dynamical times. A number n 
of phase-space positions along the orbit are recorded. Then 
the initial particle (r(to),v(to)) is replaced by n particles, 
each of which is given one of the n recorded positions and a 
statistical weight F = l/[n fs(r(to), v(to))]. Usually, out of 
20 equally sized bins in energy, we consider for splitting the 
j — 17 ... 20 bins, take n = (j — 17) 2 , and repeat the whole 
procedure a second time (with n — (j — 15)). 

(ii) fs2- ft — 0.9, ol — 1. 

This choice ensures a central density with cusp oc R~ 9 , 
which proves an effective way to increase the central number 
of particles. We do not perform any splitting. On the other 
hand the sampling is very inhomogeneous in R and z, as 
already mentioned. 

The two particles distributions arising from the two sam- 
pling distributions just described are displayed on Fig. |i[ 
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5 CUSPS IN NON-ROTATING MODELS 
5.1 Cases considered 

We assume a point mass is growing at the center of the 
model, by gas accretion, without depleting the stellar com- 
ponent (as in Young 1980, Quinlan et al. 1995, CB94, LG89). 
The BH mass grows over a time tBH from M = at t — 0, to 
the final value Mbh following (cf. Merritt & Quinlan 1998): 



M(t) 



UbhJ 

For most of the runs discussed hereafter, we took tsn = 10 7 
Yrs. 

The potential due to the BH is modeled using a Plum- 
mer potential with smoothing length tBH- 



2—) 

tBH J 



(23) 



BH 



= M/y/- 



r 2 + e 



bh- 



(24) 



We can choose esn by imposing e.g.: tBH <C Ti. For fri- 
es h < 10 _3 r/. For the standard case, where ri — 
~ 2 , we have taken £bh = 5 x 10~ 5 , and verified that 



stance, 
8 x 10" 

this value was small enough 

We consider as the standard case one where the 



5.2) 



BH 



has a mass of 2% the galaxy mass. This is supposed to be 
roughly representative of observational data. However, our 
flattened Plummer model does not have a realistic density 
profile, so that a more meaningful figure may be the ratio of 
the BH mas s to the initial core mass. The latter is roughly 
0.2 M (see §0). 



In addition, we discuss below experiments with higher 
BH masses, and different growth times. The summary of 
all the cases considered is to be found in Table ^. We will 
also consider, in the next section, models havinginitially 
some rotation. Those will be summarized in Table W of that 
section. 

We display on Fig. ^ the results for the standard case. 
The sampling distribution used was fsi- For that case, and 
other ones where the BH mass is very small, statistical noise 
is reduced by adding 10 snapshots taken within a small time 
interval after tsn- This is much less useful for the runs with 
larger BH masses (see below, Table A polar grid has been 
used to produce the graphs, with 20 points logarithmically 
spaced along r (based on Eq. ^), and 10 points equally 
spaced within [0,7r/2] along 9, which is the angle from the 
z axis. The smallest value r\ of the r- grid is chosen so that 
the total snapshot has at least 1,000 particles with r < ri. 

Quantities associated to the particles correspond to per- 
turbed quantities. The related grid quantities sli^i at grid 
points (rk,&i) are derived by applying a CIC scheme in 
the polar geometry. The total fields are finally recovered 
by adding the corresponding unperturbed analytical terms 
evaluated on the grid points. Thus for the total density: 



pki = Po(rk,9i) + piM- 



(25) 



Similarly the velocity components and dispersions are 
derived by adding the analytical terms, using Eq. |l^, for V = 
V R , V4,, and V z , and n = 1,2. Fig. § shows isolevel contours 
for the total density at the end of the simulation, as well as 
for the total velocity dispersion a 2 — a R +a 2 + a^. The dash- 
dotted curves show the initial, unperturbed quantities. Also 
displayed (first row) are the curves pt,i and Pk,w, which give 
approximately the profiles of the total density along R and 
z, respectively. The same is drawn on the bottom row for 



Table 1. Parameters for the 'standard' run 



Final M B h /Ma 

Time for BH growth 

Smoothing length for the BH 
Number of particles 



0.02 
10 7 Yrs 
5 x 10 _5 o p 
~ 300, 000 




0.5 1 1.5 2 



-2-10 1 

log(R) 



Figure 3. The two left panels show iso- contours for the den- 
sity (top) and the total velocity dispersion (bottom), for the final 
model in the standard run (solid lines). The levels are logarithmi- 
cally spaced. The initial model is shown by the dash-dotted lines 
(using different logarithmically spaced levels). The initial (dots) 
and final (full line) profiles along resp. the R axis, and the * axis, 
are also displayed, using logarithmic scales. 



the square total dispersion a 2 . The dotted segments on these 
graphs indicate a line with logarithmic slope respectively 3/2 
for the density, and 1 for the square total dispersion. 

Therefore, the figure shows that even for the less mas- 
sive BHs considered in this work, the central cusp is well 
resolved in two-dimensions. The maps also show that the 
central region has become very nearly spherical; the final 
axis-ratio is displayed on Fig. H as a function of the major- 
axis radius, scaled to rr. By eye inspection of ^, the slope 
for the velocity is very similar to what is expected theoreti- 
cally; also the density cusp appears to be similar to what is 
predicted for an initially spherical model. The cusp slopes 
will be measured with more care and discussed further in 



4 . We now turn to some runs performed in order to check 



our computations. 



5.2 Numerical checks 

A run was first made for bp — 0, in order to check that 
the slope obtained for a spherical potential agrees with the 
analytical estimates for thi s ca se. Its slope is measured, in 
the way we will explain in §5.4, both for the density and the 
velocity dispersion. We measure respectively 7 P ~ 1.52T0.02 
for the density, in good agreement with the adiabatic model 
(and the simulation by Sigurdsson et al. 1995), and 7^2 ~ 
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Figure 4. Minor to major axis ratio, as a function of R scaled 
to the influence radius rj, for the standard run (solid line). The 
dashed line shows the axis ratio for the initial model. 



0.94 ± 0.05, very close to unity, as expected in the nearly 
Keplerian potential produced around the BH. 

We also made some numerical checks for the standard 
axisymmetric run. First of all, the simulation was pursued 
after the BH growth time for an equal duration, in or- 
der to verify that a stationary model had been reached. 
The influence of the smoothing length ebh used in the 
BH potential was tested, by comparing runs using either 
a larger (cbb = 10~ 4 ap) or a smaller (cbb = 2 x 10~ 5 ap) 
smoothing length. The larger value produced a slightly shal- 
lower cusp, but results were unaffected for the smaller value. 
Therefore we use cbh = 5 x 10~ 5 ap in all subsequent runs - 
including those with larger BH masses, where a larger value 
of 6bh could have sufficed. 

We have experimented with both the fs described in 
§ [4.4 The final model obtained, for the parameters of the 
standard case, but using fs2 instead of /si, is shown on 
Fig. [B|. Results can be seen, by comparing to Fig. to be 
extremely similar to those obtained when using the other 
sampling distribution, in spite of the fact that fsi and fs2 
are very different functions. The first sampling (/si) turned 
in practice to yield results with apparently somewhat less 
particle-noise for the density. This is probably due to the 
fact that fsi corresponds to a particle distribution which 
is similar along the R and z axes, while fs2 has a cusp 
only in R. This may for instance explain why the spherically 
averaged cusp has more mass for fsi than for fs2 (see Fig. 
^|). Some perturbation mass may be less well sampled in the 
z direction. The differences, however, remain marginal. The 
two sampling function s gi ve estimates for the logarithmic 
slope of the cusp (see 5^) that can not be distinguished, 
within the error bars. Therefore we are confident that our 
results are not affected by the initial particle distribution. 

A few runs were also performed in order to make sure 
that the results do not depend significantly on numerical 
parameters, such as the spacing of the grid used for the force 
evaluation (q4.3 |), the total number of force computations by 
GRAPE (§[4.3|), or the choice that we made for the number 
of particles N = 3 x 10 5 . 
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Figure 5. Same as the previous figure, also for the parameters 
of the standard case, but using fs2 f° r the particles initialization, 
instead of fsi- The layout is as for Fig. |i[ 



5.3 Importance of the cusp self-gravity 

The perturbed central densities are high, but within small 
volumes. Therefore the mass in the cusp is not necessarily 
very large. Applying Eq. [] to the Lynden-Bell model we find 
for the influence radius: 

M BH 



ri 



2a P ~- 







(26) 



For the standard run Mbh/Mq = 0.02, thence rj 
0.16. Initially, the mass enclosed in this radius 
roughly (neglecting the effects of flattening): Mo(ri) 
Motn/Vrf + af,) 3 
lation, we find it is: 



- l (r i <r I ) Tflli 



~ 2 x 10 °M . At the end of the simu- 
M(n) = M (n) + Mi(rj) = Mo(ri) + 
2 x 10~ 5 + 6 x 10" 4 ~ 6 x 10~ 4 M . There- 
fore the mass within the cusp remains a small fraction of the 
system mass. 

On the other hand, the system has been affected by 
the central mass in a region much larger than the cusp it- 
self. This can be seen from the radial distribution of the 
perturbation mass density, displayed on Fig. ^ for the stan- 
dard case. For r > 2, dA ^ r - > < 0, corresponding to orbits 
that have been requisitioned in order to build the cusp; this 
negative perturbation extends far from the center. As a con- 
sequence, the cumulated mass Mi(r) decreases after r ~ 2. 
It is 10% of its maximum value at r ~ 5. Outside this ra- 
dius the cusp force is practically zero (Fig. ||). The orbital 
participation will be further described below (§5.5). 

We may estimate the self-consistent force (due to the 
cusp itself) F„.g. = V$i within the influence radius, where 
we assume pi(R, z) w K r -3 ^ 2 . A reasonable rough guess for 
K is such that pi(rj) ~ po(ri). This yields a roughly correct 
'calibration' for the central value, in agreement with the 
numerical value at the end of the simulation. Fig. ^| displays 
the different radial force terms: Fo — Vp&o, Fi = Vji$i, 
Fbh- Within a substantial fraction of the initial core radius 
(~ 0.75 ap — 3), the unperturbed force is negligible. The 
BH dominates the central region, within roughly ap/A, by 
a factor scaling from 10 4 at R ~ to 10 2 at R ~ rj. The 
self-consistent force contribution reaches at most ~ 1 % for 
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Figure 6. Left: The radial distribution of perturbation mass 
dMi/dr (solid line) and the cumulated distribution Mi(r) 
(dashed line), for the standard case with Mbh /Mq = 0.02, in 
percents of Mo. Right: Comparison between the forces due to the 
unperturbed mass distribution (dashed line), the BH (solid line), 
and the cusp gravity (dotted line). The insert shows the ratio of 
the cusp force F\ to the total force Fo + F\ + FbHj m percents. 




Figure 7. Same as Fig. |£j| but for a BH with mass Mbh = 0.2 Mq 

Mbh /Mo = 0.02. Therefore we have neglected it for most 
runs with this BH mass. The self-gravity contributes at most 
for ~ 2% when Mbh /Mo = 0.03 (see Table |). This is still 
very small, therefore we have also neglected the cusp self- 
gravity for models with Mbh /Mo = 0.03. 

The self-consistent term is not negligible, however, for 
the larger mass ratio Mbh /Mo = 0.2 considered in this 
work, since it contributes for up to ~ 15% of the total force 
at R~ap/A (Fig.@). 

5.4 Measure of the cusp slope 

The Poisson noise due to the discretisation into particles 
varies with the number of particles as l/y/~N. To reduce it, 
we continue the simulation for another 0.5 x 10 7 Yrs, and 
draw 10 snapshots, at regular time intervals. The snapshots 
are then merged, and this snapshot of 10 iV ~ 3 x 10 6 par- 
ticles is analyzed. This is especially useful for the runs with 
Mbh /M = 0.02, and, to a lesser extent, Mbh /Mq = 0.03. 

The central region (within ~ rr) is very nearly spheri- 
cal, as can be seen from the 2D contour map of p(R, z) (Fig. 

. Therefore, in order to measure the cusp slope, we can 
integrate over the 6 angle of the polar coordinates (r, 0). A 
logarithmic grid in r is built, using Eq. and such that the 
sphere with radius the first grid point n encloses ~ 1000 
particles. The spherically symmetrized density and velocity 
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Figure 8. Spherically symmetrized profiles for the density and 
velocity dispersions for the standard run. The total velocity dis- 
persion (thick line) and the dispersions along each direction are 
displayed. The sampling distribution used is resp. fsi (top) and 
fs2 (bottom). The positions of the filled squares indicate the grid 
used; the lines indicate slopes 7 = 3/2, for the density, and 7=1 
for cr 2 



dispersions are evaluated at the grid points, and these grid 
values used to derive the corresponding logarithmic slopes. 

The logarithmic slopes 7 P for p(r), and 7^ for cr 2 (r), 
are evaluated by computing the respective average loga- 
rithmic derivative over the grid points in some interval 
Ir = [0, r m ax]- To estimate the error in this measure, we 
also compute the slope for each individual snapshot, and 
infer the mean deviation. 

The radius r max should in theory be r max ~ rr. In 
practice, we take the maximum interval in which the radial 
profile appears as a straight line when plotted in logarith- 
mic axes. This leads actually to consider a slightly smaller 
interval than [0, rj] for c 2 (r), in order to avoid the transition 
region from cusp to core (see Fig. |^). Including grid points 
that belong to this transition region would obviously lead to 
a systematic under-estimate for the cusp slope. 

This point is illustrated in Table bl for the standard 
case with Mbh /Mo = 0.02. This is the case where the mea- 
sure depends the most on the choice of I r . Indeed, the less 
massive the BH, the fewer the grid points contained within 
rj. As a consequence, each of the grid points has a larger 
statistical weight in the average slope. 

Of course, the value of the slope for a 2 must be 1 in 
the region where the BH dominates the dynamics. As can 
be seen from the table, the effects of excluding successively 
rjv, rjv-i etc. from the N grid points within [0, n] are, for 
7 CT 2 : (i) to increase the average slope estimate (ii) to increase 
the mean deviation around the average value, as fewer data 
points are made available. 

If we exclude the 2 points lying within the cusp/core 
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Table 2. The value of the logarithmic slope, estimated for re- 
spectively p(r) and o~ 2 (r); the mean deviation is also indicated. 
The case shown is for M BH /M = 0.02, and f s = f sl . The 
successive columns correspond to taking different fractions of the 
total number N of grid points within [0, rj] —see text. 





n = N = 8 


n = N- 1 


n = N - 2 












1.46 ±0.06 


1.43 ±0.05 


1.43 ±0.06 


7 CT 2 


0.93 ± 0.05 


0.99 ±0.06 


1.01 ±0.08 



transition, the logarithmic slope measured for a is very 
nearly equal to unity, as it should be. 

As for the density, there are up to 8 points that lie well 
within the cusp (Fig. |8|). Our density profile deviates slightly 
from a pure power-law at the 2 inner points. Therefore we 
obtain a better estimate of the slope in the linear section 
by taking into account the full set of points within rj. With 
these grid points, we find for the standard case: 

• using /si : 7 P ~ 1.46 ± 0.06 

Fig. ^ also displays the density and velocity profiles 
obtained when the sampling distribution is fs2- The cen- 
tral density value is slightly smaller than with the previous 
case, probably because the sampling is not as efficient due 
to its strong spherical asymmetry. The value measured for 
the density logarithmic slope is now: 

• using f S 2 ■ 7p — 1-35 ± 0.15, 

This value is somewhat smaller than the value for the 
other sampling, but still consistent with it. We therefore con- 
clude that the two sampling distributions tested, in spite of 
their different behaviour in the centre, yield similar results. 
We therefore believe the results are not an artifact of the 
method used. The value of y p is consistent, within its error 
bars, with the value given in the literature for the spherical 
case, 7 = 3/2. We find no evidence for a significantly larger, 
nor smaller, slope. 

The results found for exp eriments using a larger final 
BH mass are discussed in §5.7. For such cases, the dispersion 

as the number 



in the measures is very small (see Table 
of points useful for the measure is larger. 



5.5 Orbital response to the BH 

When we use the PP method, we have direct access to the 
distribution functions fo, fi and f(E, L z ) = fo + f\. Those 
are available with a large signal-to-noise ratio, even for the 
case where the BH mass is only 1% of the model mass, that 
we analyse in this paragraph. 

The initial model can be viewed by plotting the surface 
fo(Eo, Lz), with Eq = $o — v 2 /2 the unperturbed energy. 
Since $o(r = 0) = 0.25, < E < 0.25. This surface, drawn 
from a set of initial conditions for the particles, is shown on 
Fig. |^ (central panel). We use a regular grid in Eq and L z to 
bin the particles, and sum the masses ma = fo(wi) / fs(u>i) 
within each cell of the grid. It could of course have been 
drawn analytically, but, as it is, the figure shows how smooth 
the numerical realization is. A map using a linear grey scale 
is shown on top of the surface. At fixed energy, the d.f. in- 
creases with increasing L z , as can be seen from the figure. 



This is related to the excess of tangential motion supporting 
the flattening. 

In a similar way, using snapshots at the end of the stan- 
dard run, we can build the surfaces fi(E, L z ) and f{E, Lz). 
We plot in fact two slightly different maps, (i) A map 
of /i as a function of the initial values of Eo and L z - 
L z is anyway conserved during the BH growth- reveals 
which orbits have been most perturbed, (ii) On the other 
hand, a map for f(E,L z ) — fo(E,L z ) + fi(E,L z ) (with 

E = $o + <&bh + $s. g . — \ the total final energy) yields 
information about the orbital structure in the final model. 
In practise and for the sake of simplicity , w e neglect the 
self-gravity of the cusp (we have shown in 



5.3 that the self- 



gravity contribution is small for Mbh/Mo < 0.03), and plot 
f(E*,L z ), where E* = $ + $>bh - 

The map for f\ is show on Fig. ^| (left panel). Pos- 
itive values of fi are displayed in grey shades darker than 
the background (with white contour lines), whereas negative 
values appear as lighter shades (with black contour lines). 
The map shows that orbits having Eo ~ $o have been the 
most perturbed. Also, at a given energy, the perturbation 
is more important for orbits with small angular momentum 
\L Z \. For the most bound orbits (Eo — $o), fi > 0, reflect- 
ing the fact that they now belong to phase-space regions 
more populated than initially. By contrast, a negative per- 
turbation (f\ < 0) can be found for smaller values of Eo, 
and L z ~ 0. These orbits have an apocenter located at a 
radius much larger than the influence radius. The fact that 
the perturbation extends much fu rthe r than the influence 
radius was already seen in Fig. ^ (§5.3). 

We have used 10 snapshots, issued at the end of the sim- 
ulation and merged together, in order to produce the surface 
f(E*, L z ) shown on the right panel of Fig. ^j. The final total 
distribution shows most clearly that only orbits with L z ~ 
have a final energy E* > <f>o(r = 0) (Fig. W. These are orbits 

2 

such that &bh — \ > 0, i.e. orbits that have become bound 
to the BH potential. Such orbits must provide an important 
contribution to the cusp. However, the fact that nearly ra- 
dial orbits are more efficiently attracted into the cusp, does 
not imply radially biased motion within the cusp, since these 
orbits get rounder (see also CB94). In models with a cen- 
tral cusp, constructed either with, or without, a BH, Dehnen 
(1995) and Dehnen & Gerhard (1994) find that such orbits 
typically have L z ~ L C (E). 

The final d.f. appears remarkably constant at high bind- 
ing energies (f(E) ~ est for E > <l?o(0) and Lz ~ 0). 
Tremaine et al. (1994) have built analytical models for spher- 
ical systems with a central density cusp and a central BH. 
For a cusp with slope 7 = 3/2, their d.f. tends to a con- 
stant for high binding energies, in a way very similar to our 
results. 

The final d.f. can therefore be said 'degenerate' at high 
binding energies, in the sense that different energy levels 
have the same probability. Such a degeneracy is also typical 
of violently relaxed systems (see Lynden-Bell 1967, Chava- 
nis & Sommeira 1998). Our results therefore support those 
by Stiavelli (1998), who has shown that the observable sig- 
natures of an adiabatically grown BH, on the one hand, and 
of violent relaxation around a pre-existing BH, on the other 
hand, are very similar. We have shown that this is due to the 
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Table 3. List of experiments with non-rotating models, with their numerical results. Column 4 gives the logarithmic 
slope measured for the density, and column 5 gives that for the total velocity dispersion. 1 1 Mi 1 1 gives an estimate of 
the mass in the perturbation. F max is the maximum relative contribution to the total force of the self-gravitv force. 
The last column gives an estimate of the relative error made in evaluating the perturbed quantities (see §4.1). 
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Figure 9. Distribution functions for the case without net rotation. The left panel shows the perturbation distri- 
bution /i(i?o, L z ) as a function of the initial unperturbed energy and angular momentum of the particles. The grey 
scale is dark for positive values, scaling linearly from 10 — 5 to 4x 10 — 3 (with white contours over-imposed), and light 
(with black contours) for negative values. The 3-D plots display the initial fo(Eo, L z ) (center) and final f(E,L z ) 
(right) distribution functions- see text for details on the corresponding E, L z axes. 



fact the distribution functions themselves have very similar 
behaviour at high energies. 



5.6 Influence of growth time 

The influence of the BH growth time on the final cusp has 
been checked, by experimenting with shorter times (Tbh = 
10 6 , 10 5 Yrs), for the case M BH /M = 0.03. The dynamical 
times within the cusp are close to 10 Yrs, therefore the 
assumption of adiabaticity does not hold for the latter case. 
Moreover, some orbits may not have had time to rearrange 
into an equilibrium during the BH growth. Therefore all 
simulations are pursued after the BH has reached its final 
mass, so that the total time is 10 7 Yrs for all cases. 

Results for the logarithmic slopes are reported in Table 
^| In good agreement with Sigurdsson et al. (1995), we found 
that the adiabatic model is still roughly valid for growth 
times of order ~ 10 Td yn - Nevertheless small growth times 
introduce some differences, and already for Tbh = 10 6 Yrs, 
there is a tendency for the cusp to spread over a larger di- 
ameter. Whereas the density profile follows a power-law in a 
smaller region around the center, the transition zone around 
77 is much larger than in the case Tbh = 10 7 Yrs. There is 



therefore a trend for the average cusp slope to be slightly 
smaller than 3/2. 

This trend towards a shallower cusp is more evident 
in the case where Tbh = 10 5 Yrs (Fig. [u]). Actions are 
not expected to be conserved in this case where T grow ~ 
Tdy n , as orbits are deflected by the rapidly varying central 
potential. Probably as a result of such deflections, we find 
an excess of radial velocity dispersion (or ~ o z > a^) in the 
central region, as found also by Sigurdsson et al. (1995) for 
experiments with a similar growth time. Less orbits settle in 
the vicinity of the BH, and the cusp appears less strong than 
in the adiabatic case, with a logarithmic slope 7 P ~ 0.9. 

5.7 Varying the BH mass 

We also experiment with BH masses which are a larger frac- 
tion of the model mass (see Table g). The mass of the dark 
matter concentrations betrayed by a central cusp in early- 
type galaxies is evaluated to be at most 2% of the galaxy 
mass (e.g. Richstone et al. 1998). Maybe more significant for 
the cusp morphology and kinematics, is the ratio between 
the BH mass and the initial core mass of the galaxy. 

Semi-analytic works (Young 1980, Quinlan et al. 1995) 
predict that the density profile of the cusp remains roughly 
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Figure 10. Density profile for a model with final Mbh /Mq = 
0.03 and with BH growth time respectively Tbh = 10 6 Yrs (left) 
and Tbh = 10 5 Yrs (right). The line segments indicate the 7 = 
3/2 power-law. 
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Figure 12. Left: Isocontours for the final density, and profiles 
along R and z, for Mbh /Mo = 0.2. Isocontours and profiles for 
the velocity dispersion are displayed on the second line. 
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Figure 11. Spherically symmetrized profiles for p(r) and <r 2 (r) = 
f/3(cr|j + a 2 z + ct?), shown for Mbh /Mcore = 0.15 (top panels) 
and Mbh /Mcore = 1 (bottom panels). The dispersions along 
individual directions a?., a\, and a 



R 



7 ^ have a shape very similar to 



that of the total velocity dispersion. 



self-similar, as long as the final values of the BH mass are 
smaller than the mass within the initial core. On the other 
hand, for Mbh larger than the core mass, the cusp is steeper, 
with 7 P approaching a value of 2 when the mass ratio is 
Mbh / M core 1. Of course no variation is seen in the ve- 
locity dispersion cusp, which is always what we expect in an 
essentially keplerian potential (7 CT 2 = 1). 

The cases reported in Table |3| correspond to a BH with 
mass equal respectively to 0.05, 0.15 and 1 times the core 
mass M cor e- The first case has already been discussed in 
§5.4; the two other cases are illustrated on Fig. [ll], which 
displays the spherically symmetrized profiles for the density 
and the velocity dispersions. 

As can be seen from this Figure, the number of points 



within the power-law cusp is much larger than for the stan- 
dard case. For a BH whose mass is a substantial fraction of 
the core mass, the statistics on the results is therefore much 
better than in the standard case Mbh / M CO re = 0.1. 

A single snapshot (with N ~ 320, 000 particles) suffices 
to trace clearly the cusp, and measure its slope with good 
accuracy. The dispersion in this measure is very small, as can 
be seen in the results summarized in Table ^. Fig. [l2| shows 
in more detail the case with Mbh = 0.2 Mo, corresponding 
to Mbh / Mcore — 1. It was plotted using one snapshot. 

The logarithmic slope 7 P that we measure for the total 
density is an increasing function of the final BH mass (see 
Table |^). It takes values within [3/2,2], very much as ex- 
pected from the spherical adiabatic model. Fig. |l^ compares 
the density profiles obtained for the different Mbh /M CO re 
values considered. 

We also ran a simulation where the final BH mass was 
2 Mcore- Although results seemed very reasonable in all the 
cases considered, the error in the mass conservation increases 
in this case to non-negligible values (Mi / j | A/i j | ~ 20%) . The 
computation is therefore at the limit of credibility for the 
method. The perturbation technique used here is indeed bet- 
ter suited to small or moderate values of Mbh /Mo, rather 
than larger ones. The reason is not that the method is lim- 
ited to the linear regime: the computation is fully non-linear, 
since we compute perturbed orbits within the full poten- 
tial. The reason is rather that the sampling is more difficult 
to control for larger perturbations. The non-conservation of 
the total mass reflects the sampling inadequacies. For higher 
mass ratios than those considered in this work, a standard 
iV-body technique should probably be preferred. 



6 CUSPS IN ROTATING SYSTEMS 

In this section, we investigate the effects of an initial rotation 
of the host galaxy. LG89 and CB94 both study the rotation 
brought to central regions by a growing cusp, but find dif- 
ferent results. Is the galaxy response somehow affected in 
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Figure 13. Density profile along r obtained for a BH with final 
mass resp. 0.05,0.15, 1. times the initial core mass 



Figure 14. Initial rotation velocity profile along the Ft axis, for 
the model with resp. p = pi (solid line), and p = p2 (dashed line), 
as given by Eq. 



these works by the approximation of a spherical potential? 
Here we can study this case in its true geometry, without 
any such approximation. Moreover, we would like to clarify 
which difference in these works is responsible for the differ- 
ence in the final velocity field. The net rotation within the 
influence radius is weak in CB94, unless the BH has a mass 
much larger than the initial core mass. On the other hand, 
LG89 obtain a cusp in the rotation velocity for any initial 
BH mass. We therefore investigate in turn models that are 
built similarly to those considered by CB94, and LG89. 

6.1 Models with pi(L z ) 

We first experiment with rotating models built using 
p{E,L z ) = pi(L z ) (Eq. pd| . The motivation for consider- 
ing this kind of rotation tapering for small L z , is obviously 
its simplicity. Furthermore, such models are not unplausible 
if a galaxy has formed through a major merger. Violent re- 
laxation in the central regions can erase net rotation for very 
bound orbits, whereas less bound orbits may still exhibit a 
preference for one sense of rotation. An example of merger 
remnant with this sort of angular momentum distribution 
is found in numerical simulations by Barnes & Hernquist 
(1996; their Fig. 17). 

The parameters in pi have been chosen as follows: 

L m = 0.2Z V = 0.5. (27) 

The corresponding rotation curve is displayed on Fig. [l4|. 

For each BH mass considered, the slope of the cusp is 
the same, within the error bars, as in the corresponding non- 
rotating case. These results are summarized in Table m 

The initial distribution fo(Eo,L z ) is displayed on the 
middle panel of Fig. |l5| for a mass ratio Mbh /Mo = 0.01, 
corresponding to Mbh /M core ~ 0.05. Again, even for this 
small final BH mass, the distribution function exhibits little 
noise. Preference for the positive sense of rotation is reflected 
by the distortion of the surface fo(Eo, L z ), which is higher 
for L z > values. 

The left panel of Fig. [l^ shows the perturbation /i as a 
function of the initial integrals of motion Eo,L z . The grey 
shades and contour lines are as in the left panel of Fig. bl 



Again, perturbation is strongest for large Eo values, and 
small L z values. The peak of positive perturbation (at E ~ 
Eo) is roughly even in L z , corresponding to an energy range 
where the initial distribution fo is itself roughly even in L z . 
We have indeed considered a model like the model in CB94, 
where the net rotation is mainly due to orbits having L z > 
Lo, with Lo some finite value. For orbits with high binding 
energy, L C (E) is smaller than Lo, so that practically none 
of them can contribute to the net rotation. On the other 
hand, for smaller energies Eq, the perturbation fi is odd 
with respect to L z . This is a consequence of the odd term 
rip(L z ) that has been added to the distribution function fo 
in order to construct fo - see Eq. []. 

The final total distribution function f(E,L z ) is dis- 
played on the same figure (right panel). At high energies, 
the distortion of this surface towards positive L z has slightly 
decreased in amplitude, corresponding to the negative per- 
turbation /i for the positive L z . The strip of particles having 
E > $o(0) is very narrow around L z = 0, as can be best 
seen on the grey scale map on top of the right panel. Thus, 
only orbits with L z ~ become bound to the cusp. As a 
consequence, orbits confined to the central regions do not 
create a significant global rotation. 

Fig. hq shows the rotation at the end of the BH growth 
for different mass ratios. Similar to results by CB94, we find 
that a moderate BH mass does not bring a significant ro- 
tation within the influence radius. This result therefore was 
not affected by simplifications in the analytical derivation of 
CB94. Our results are indeed very close to those of CB94, 
as can be seen by comparing our Fig. |l^ and their Fig. 3. 

As the BH mass increases well above the core mass, its 
influence radius eventually overcomes the radius correspond- 
ing to the circular orbit with angular momentum Lo, and 
rotation becomes more important in the centre. However, 
for the mass ratios we consider (up to Mbh / M core ~ 1), we 
find that no cusp is produced in the rotation velocity. 

6.2 Models with p 2 {E,L z ) 

We now experiment with a rotating model that uses 
p{E,L z ) = p 2 {L z /L c {E)) given by Eq. [TTJ. This model is 
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Table 4. List of numerical experiments with rotating models. See Table |3j for the symbols definition. The results 
are shown for p = pi(E, L z )\ experiments with p = P2(E, L z ) yield practically identical measures for the logarithmic 
slopes of both p and a 2 . 
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Figure 15. Distribution functions for the case Mbh /Mq = 0.01 with rotation, and p = p\. The three panels are as in 
Fig. ^, displaying respectively f±(Eo, L z ) (left), faiEo, L z ) (middle) and f(E, L z )(right). The grey scales and contour 
lines used for the left panel are defined as in Fig. H 
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Figure 16. Rotation velocity profiles, along the R axis, for three 
different mass ratios Mbh /Mcore = 0.05,0.15, 1, as labeled; the 
dotted line shows the initial rotation of the model 



similar to that by LG89, with their parameter uj correspond- 
ing roughly to n/(2x) in our models. We obtain a model with 
an initial rotation very close to the rotation for the models 
with pi previously considered, by taking: 



0.2, 



0.5. 



(28) 



The corresponding rotation profile along the R axis can be 
seen on Fig. 

In the models now considered, an identical fraction of 
orbits, at any energy level, contributes to the net rotation. 




0.05 0.10 0.15 0.20 0.25 
E 



Figure 17. Perturbation distribution function /i as a function 
of the initial energy Eq and L z , for the case Mbh /Mq = 0.01 
with rotation, and p = pi. 



Orbits having initially E « 3>o(0) have L C (E) ss 0, so that 
even orbits with very small values of L z contribute to the 
net rotation. These orbits are those which will bring rot atio n 
A map for fi(Eo(t = 0),L Z ) (see 



within the cusp. A map for fi(Eo(t = 0),L Z ) (see §5.5), 
displayed on Fig. shows indeed that, for orbits with high 
initial energies, the perturbation has grown positive for L z > 
0, and negative for L z < 0. Therefore a rotational velocity 
has efficiently built within the cusp. 

The effects on the final distribution function are best 
seen from Fig. [ll]. We have plotted the case withMs// /Mo — 
0.03, rather than the standard case with Mbh /Mo = 0.01, 
because the structure within the total final d.f. is more easily 



© 2000 RAS, MNRAS 000, §4l^ 



15 




ies (see Jaffe et al. 1994, and Lauer et al. 1995). Such a 
disk would of course demand a different model for the ro- 
tation profile. The suggestion that the power-law galaxies 
(those with steep cusps) owed their inner luminosity profile 
to the presence of a disk has been much weakened by the 
fact that Lauer et al. (1995) find little evidence for inner 
disks in the power-law E's of their sample. More evidence 
for central disks, at the scale of ~ 10 pc, has been collected 
for SO galaxies (van den Bosch et al. 1994, Scorza & van den 
Bosch 1998), or giant E's with shallow cusps (Forbes 1994). 
Obviously the statistics on the detection of such disks is still 
too small, and the dynamical models far from entering such 
detail. 



Figure 18. Final total distribution functions for the case with 
rotation, and p = p2\ the case with Mgjj / Mq = 0.03 is shown, 
for greater legibility. 




Figure 19. Same as Fig. but for an initial model with 
p(E,L z ) = p2 - see text. The label next to each curve refers 
to the mass ratio Mbh /M core . 

visible. The fact that rotation has been dragged within the 
density cusp is disclosed by the bent of the surface strip of 
particles with high binding energies. 

This explains why a cusp in rotation velocity is now 
produced for any BH mass considered, in agreement with 
the models described by LG89. The cusp follows from the 
fact that, initially, rotation was present at all levels of energy, 
including high energies. The rotation velocity profiles along 
the R axis are shown on Fig. JH| for different masses of the 
BH. The different curves are labeled as in Fig. |l6| 

Therefore we have shown that the rotation in the cusp 
region depends on the orbital structure that existed in the 
galaxy prior to the BH growth. The formation of a cusp in 
Vrot around a BH that grows adiabatically is sensitive to 
the way the initial d.f. depends on L z . If rotation is not 
present initially at high binding energies, only very massive 
BHs -with respect to the initial core mass - will produce an 
observable signature onto the V ro t profile. 

There has been some debate about the role played by a 
nuclear disk onto the light profile in the center of E galax- 



7 CONCLUSION 

This paper aims at investigating the cusp induced in ellipti- 
cal galaxies by the growth of a central, supermassive black 
hole. The spherical case is the only one to have been stud- 
ied in detail in the literature previously, using an adiabatic 
model. In this paper we study the cusp due to a growing 
BH by numerical means, in order to free ourselves from the 
assumption of spherical symmetry. 

We thus investigate the possible influence of flattening 
and rotation, by considering a simple axisymmetric model. 
For our models, which have a substantial central flattening, 
we do not find any sizeable signature on the cusp slope with 
respect to the spherical case. We deduce from this result that 
a reasonable degree of tangential anisotropy in the stellar ve- 
locities, which sustains the flattening, has little effect on the 
cusp slope. The spherical adiabatic models appear very ro- 
bust with respect to the geometry of the initial galaxy. We 
find a cusp slope of 1.5 in the models where the BH mass 
is less than the initial core mass, whereas the observed pre- 
ferred value for low mass E's is larger, around 1.9. This fact 
remains to be explained. One possibility is that the BH grew 
in a galaxy with initially a very small core, as suggested by 
van der Marel (1999), or even no core at all. The obser- 
vational trends (steeper cusps for low luminosity E galax- 
ies, and smaller cores for less massive core Ellipticals) are 
roughly consistent with this suggestion (see van der Marel 
1999). Violent relaxation during a merger or a gravitation- 
nal collapse favours highly concentrated systems (Farouki 
et al. 1983; see also Chavanis & Sommeira 1998), specially 
if dissipation occurs (Udry 1993). Gas or stellar collisions 
around the BH may also have favoured steeper cusps in low 
L galaxies. The relaxation time is indeed shorter for less 
massive E galaxies as they have denser centers (cf. Magor- 
rian & Tremaine 1999) and we know that the cusp induced in 
a collisional system is higher than in the collisionless models, 
with 7 P around 7/4 (cf. Duncan & Shapiro 1983). 

The formation of a cusp in V ro t around a BH that grows 
adiabatically appears sensitive to the way the initial d.f. 
depends on L z . If rotation is not present initially at high 
binding energies, only very massive BHs will produce a sig- 
nature in the V ro t profile. Such an initial orbital distribu- 
tion could for instance be expected after mergers involving 
an efficient violent relaxation in the central region. Stellar 
kinematic observations with the required resolution are ex- 
tremely scarce, although STIS observations {e.g. Joseph et 
al. 2000) may soon improve on this situation. Existing data 
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for M32 (Joseph et al. 2000) and for the E5 galaxy NGC3377 
(Kormendy et al. 1998) show an increase of rotation velocity 
within the inner few arcsec. This is again consistent with the 
interpretation that in power-law elliptical galaxies the BH 
mass exceeds the initial core mass (eventually zero). 

On the other hand, a cusp in rotation is produced, what- 
ever the BH/core mass ratio, when initially the fraction of 
orbits contributing to rotation is non-zero at high binding 
energies. Evidence for such a cusp is probably present in 
HST observations of the SO galaxy NGC 4342 (Cretton & 
van de Bosch 1999). 

Let us finally note that the spherical adiabatic model 
yields values for the mass of dark matter concentration sim- 
ilar to more elaborate dynamical modeling (van der Marel 
1999). This however does not preclude that central BHs pre- 
existed to their host galaxy formation, since the adiabatic 
growth scenario predicts final galaxy models very similar to 
those produced by violent relaxation around a pre-existing 
BH (as demonstrated by Stiavelli 1998, and by our own re- 
sults on the distribution function). 
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